      SUBROUTINE  RPRSWF(C,XI,M,N,MINT,INDX,NMAX,MAXJW,JW,D,ALAMD,      00010000
     +                   SIGMA,DSUM,R1,DR1,R2,DR2)                      00020012
CC                                                                      00030000
CC        THE PROLATE SPHEROIDAL RADIAL FUNCTIONS WITH COMPLEX ARGUMENTS00040000
CC        THE SIZE PARAMETER C IS REAL   -   REAL FUNCTIONS             00050000
CC        INDX=1,  EXTERNAL FIELD*C1   ;   INDX=2,  INTERNAL FIELD*C2   00060000
CC        MINT = INDEX PARAMETER FOR CALCULATION OF BESSEL FUNCTIONS.   00070000
CC        M .LE. 40     :     JW .LE. 150                               00080011
CC        THE SECOND KIND FUNCTIONS ARE NEVER CALCULATED FOR INDX = 2.  00090000
CC     -----------------------------  REVISED IN 1978 ----------------- 00100000
      IMPLICIT REAL*8(A-H,O-Z)                                          00110000
      PARAMETER  (JMX=150,MXPQ=190,MXSJ=150,MXSN=190)                   00120011
      DIMENSION  DNRAT(80),DN(80),DPRAT(200),DP(200),AMX(JMX),D(JMX),   00130011
     1           SJ(MXSJ,2),DSJ(MXSJ,2),SN(MXSN,2),DSN(MXSN,2),         00140003
     2           PX(MXPQ),DPX(MXPQ),QX(MXPQ),DQX(MXPQ),QN(80)           00150011
      FLOAT(N)=DFLOAT(N)                                                00160000
      E(J)=C2*DFLOAT((M+M+J)*(M+M+J-1))/DFLOAT((2*(M+J)-1)*(2*(M+J)+1)) 00170000
      F(J)=DFLOAT((M+J)*(M+J+1))+C2*DFLOAT(2*(M+J)*(M+J+1)-2*M*M-1)/    00180000
     +     DFLOAT((2*(M+J)-1)*(2*(M+J)+3))-ALAMD                        00190003
      G(J)=C2*DFLOAT((J+1)*(J+2))/DFLOAT((2*(M+J)+1)*(2*(M+J)+3))       00200000
C                                                                       00210000
      XI2=XI*XI-1.D0                                                    00220000
      C2=C*C                                                            00230000
      H= C                                                              00240000
      HXX=H*XI*XI                                                       00250000
      LMTN=0.835*HXX                                                    00260000
C                                                                       00270014
        IF(M.EQ.0)  THEN                                                00280003
      XIA=1.D0                                                          00290003
      XIB=0.D0                                                          00300003
        ELSE                                                            00310003
      XIA=DSQRT(XI2/(XI*XI))**M                                         00320000
      XIB=DFLOAT(M)/(XI*XI2)                                            00330000
        END IF                                                          00340003
      L=N-M                                                             00350003
      PSCRT=(1.0E-12)*(10.0**(L/10.0))                                  00360003
       IF(MAXJW.GE.JMX-1)  MAXJW= JMX-2                                 00370011
      R2=0.D0                                                           00380000
      DR2=0.D0                                                          00390000
C                                                                       00400003
C         SPHERICAL BESSEL FUNCTIONS                                    00410000
C                                                                       00420003
      IF(M.EQ.MINT.AND.N.EQ.M) GO TO 4002                               00430000
      GO TO 4005                                                        00440000
 4002 NSJ=MAXJW+M                                                       00450000
       IF(NSJ.GE.MXSJ-1)  NSJ= MXSJ-2                                   00460013
      DX=C*XI                                                           00470000
      NUP=10+DX/5.                                                      00480000
      NSWT=NSJ                                                          00490000
      NA=NSJ+NUP                                                        00500000
      W=( DSIN(DX)-DX* DCOS(DX))/(DX*DX)                                00510000
      T3= 0.D0                                                          00520000
      T2= 1.D-70                                                        00530000
      DO 4000 NB=1,NA                                                   00540000
      NC=NA+1-NB                                                        00550000
      T1=DFLOAT(NC+NC+1)*T2/DX-T3                                       00560000
      IF( DABS(T1)-1.D+70) 4001,4003,4003                               00570003
 4003 T1=T1*1.0D-70                                                     00580000
      T2=T2*1.0D-70                                                     00590000
      NSWT=NC                                                           00600000
 4001 T3=T2                                                             00610000
      T2=T1                                                             00620000
      IF(NC.GT.NSJ) GO TO 4000                                          00630000
      SJ(NC,INDX)=T1                                                    00640000
 4000 CONTINUE                                                          00650000
      RATIO=W/SJ(2,INDX)                                                00660000
      DO 4010 ND=1,NSJ                                                  00670000
      IF(ND.GT.NSWT) GO TO 4015                                         00680000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)                                     00690000
      GO TO 4010                                                        00700000
 4015 SJFACT=1.D-70                                                     00710000
      IF( DABS(SJ(ND,INDX)).LT.1.0D-5) SJFACT=0.D0                      00720000
      SJ(ND,INDX)=RATIO*SJ(ND,INDX)*SJFACT                              00730000
 4010 CONTINUE                                                          00740000
      DSJ(1,INDX)=-SJ(2,INDX)                                           00750000
      DO 4020 NE=2,NSJ                                                  00760000
      NF=NE-1                                                           00770000
 4020 DSJ(NF+1,INDX)=SJ(NF,INDX)-DFLOAT(NF+1)*SJ(NF+1,INDX)/DX          00780000
C                                                                       00790003
C         SPHERICAL NEUMANN FUNCTIONS                                   00800003
C                                                                       00810003
      IF(INDX.EQ.2) GO TO 4005                                          00820000
      IF(XI.LT.1.005) GO TO 4005                                        00830000
      IF(HXX.LT.10.0) GO TO 4005                                        00840000
      SCLSN=1.0E-20                                                     00850000
      IF(DABS(DX).GE.15.0) SCLSN=1.0E-10                                00860000
      SN(1,INDX)=- DCOS(DX)/DX                                          00870000
      SN(2,INDX)=-( DCOS(DX)+DX* DSIN(DX))/(DX*DX)                      00880000
      SN(1,INDX)=SN(1,INDX)*SCLSN                                       00890000
      SN(2,INDX)=SN(2,INDX)*SCLSN                                       00900000
      DSN(1,INDX)=-SN(2,INDX)                                           00910000
      MAXSN=MAXJW+M                                                     00920000
       IF(MAXSN.GE.MXSN-1)  MAXSN= MXSN-2                               00930011
      DO 4030 NP=2,MAXSN                                                00940000
      NQ=NP-1                                                           00950000
      IF(NQ.LE.1) GO TO 4040                                            00960000
      SN(NQ+1,INDX)=DFLOAT(NQ+NQ-1)*SN(NQ,INDX)/DX-SN(NQ-1,INDX)        00970000
 4040 DSN(NQ+1,INDX)=SN(NQ,INDX)-DFLOAT(NQ+1)*SN(NQ+1,INDX)/DX          00980000
 4030 CONTINUE                                                          00990000
 4005 CONTINUE                                                          01000000
C                                                                       01010003
C         ASSOCIATED LEGNDRE FUNCTIONS                                  01020000
C                                                                       01030003
      IF(INDX.EQ.2) GO TO 4007                                          01040000
      IF(N.GT.M) GO TO 4007                                             01050000
      MAXPQ=MAXJW+M                                                     01060000
       IF(MAXPQ.GE.MXPQ-1)  MAXPQ= MXPQ-2                               01070011
         CALL ALEGQ(XI,M,MAXPQ,QX)                                      01080003
      MLEG=0                                                            01090000
      ALEG=1.D0                                                         01100000
 4070 MLEG=MLEG+1                                                       01110000
      ALEG=ALEG*DFLOAT(2*MLEG-1)                                        01120000
      IF(MLEG.LT.M)  GO TO 4070                                         01130000
      PX(1)=ALEG*(DSQRT(XI2)**M)                                        01140000
      PX(2)=DFLOAT(2*M+1)*XI*PX(1)                                      01150000
      DO 4050 I=1,MAXPQ                                                 01160000
      MI=M+I                                                            01170000
      PX(I+2)=(DFLOAT(MI+MI+1)*XI*PX(I+1)-DFLOAT(M+MI)*PX(I))/FLOAT(I+1)01180000
      DPX(I)=(DFLOAT(I)*PX(I+1)-DFLOAT(MI)*XI*PX(I))/XI2                01190000
      IF(I.EQ.1) GO TO 4050                                             01200000
      DQX(I-1)=(DFLOAT(I-1-M)*QX(I)-DFLOAT(I-1)*XI*QX(I-1))/XI2         01210000
 4050 CONTINUE                                                          01220000
      IF(M.EQ.0) GO TO 4007                                             01230000
      S1=QX(2)                                                          01240000
      S2=QX(1)                                                          01250000
      DO 4080 KA=1,M                                                    01260000
      KK=KA-1                                                           01270000
      S3=(DFLOAT(1-2*KK)*XI*S2+DFLOAT(KK+M-1)*S1)/DFLOAT(M-KK)          01280000
      S1=S2                                                             01290000
      S2=S3                                                             01300000
      QN(KA)=S3                                                         01310000
 4080 CONTINUE                                                          01320000
 4007 CONTINUE                                                          01330000
C                                                                       01340003
C         EIGENVALUES AND EXPANSION COEFFICIENTS  -  D                  01350000
C                                                                       01360003
      IF(N.GT.M) GO TO 4009                                             01370000
         CALL DMATR(M,C,NMAX,AMX,1)                                     01380003
 4009 MTRX=L+1                                                          01390000
      AINT=AMX(MTRX)                                                    01400000
C                                                                       01410013
          CALL  RDCOEF(H,M,N,1,AINT,JW,D,ALAMD,PI,SIGMA,DSUM,FCTR)      01420012
C                                                                       01420112
      IF( DABS(DSUM).LT.1.0D-20) GO TO 3                                01430006
          IF(INDX.EQ.1)  THEN                                           01440015
C*       WRITE(6,777)  H,M,N,JW,AINT,ALAMD,PI,SIGMA,DSUM,FCTR           01440115
C*777    FORMAT(1H , F7.3, 3I5,1P6E15.4 )                               01441015
          END IF                                                        01442015
C                                                                       01443015
CC            RADIAL FUNCTION AND ITS DERIVATIVE - RJ(M,N,HX),DRJ       01450003
CC            THE RADIAL FUNCTIONS OF FIRST KIND IN SPHERICAL BESSEL    01460003
      IFIN=JW                                                           01470000
      IF(MOD(L,2).EQ.0) IFIN=JW+1                                       01480000
      R1A= 0.D0                                                         01490000
      DR1A= 0.D0                                                        01500000
      DO 1060 JR=1,IFIN,2                                               01510000
      J=JR-1                                                            01520000
      IF(MOD(L,2).EQ.1) J=JR                                            01530000
      MR=M+J                                                            01540000
      IF(MR.GT.NSJ) GO TO 1060                                          01550000
      IR=IABS(J-L)                                                      01560000
      Q=-1.D0                                                           01570000
      IF(MOD(IR,4).EQ.0) Q=1.D0                                         01580000
      RFAC1=1.D0                                                        01590016
      RFAC2=1.D0                                                        01600013
       IF(M.GE.1)  THEN                                                 01610013
      DO 15 IP=1, M                                                     01620013
      RFAC1= RFAC1*DFLOAT(J+IP)                                         01630013
   15 RFAC2= RFAC2*DFLOAT(J+IP+M)                                       01640013
       END IF                                                           01650013
   13 RCONV1= D(JR)*RFAC1*DSJ(MR+1,INDX)*RFAC2/FCTR                     01660016
      R1A=R1A+Q*D(JR)*RFAC1*SJ(MR+1,INDX)*RFAC2/FCTR                    01670016
      DR1A=DR1A+Q*RCONV1                                                01680000
      IF( DABS(RCONV1/DR1A).LE.1.0D-15) GO TO 5                         01690000
 1060 CONTINUE                                                          01700000
    5 R1= XIA*R1A/PI                                                    01710016
      DR1= XIB*R1+C*XIA*DR1A/PI                                         01720016
      IF(INDX.EQ.2) GO TO 4                                             01730000
C                                                                       01740003
C             RADIAL FUNCTION OF SECOND KIND IN TERMS OF SPHERICAL NEUMA01750000
C             SERIES WITH INTEGRAL METHOD BY SINHA & MACPHIE(1974)      01760000
      DFN=1.D+10                                                        01770000
      R2N= 0.D0                                                         01780000
      DR2N= 0.D0                                                        01790000
      IF(XI.LT.1.005) GO TO 7                                           01800000
      IF(HXX.LT.10.0) GO TO 7                                           01810000
      IF(L.GT.LMTN) GO TO 7                                             01820000
      IFN=IFIN                                                          01830000
    6 IFN=IFN-2                                                         01840000
      NR=IFN-2                                                          01850000
      R2A=0.D0                                                          01860000
      DR2A= 0.D0                                                        01870000
      DO 1070 JR=1,IFN,2                                                01880000
      J=JR-1                                                            01890000
      IF(MOD(L,2).EQ.1) J=JR                                            01900000
      MR=M+J                                                            01910000
      IF(MR.GE.MAXSN) GO TO 6                                           01920000
      IR=IABS(J-L)                                                      01930000
      Q=-1.D0                                                           01940000
      IF(MOD(IR,4).EQ.0) Q=1.D0                                         01950000
      RFAC3= 1.D0                                                       01960013
      RFAC4= 1.D0                                                       01970013
        IF(M.GE.1)  THEN                                                01980013
      DO 17 IP=1, M                                                     01990013
      RFAC3= RFAC3*DFLOAT(J+IP)                                         02000013
   17 RFAC4= RFAC4*DFLOAT(J+IP+M)                                       02010013
        END IF                                                          02020013
      RCONV2= Q*D(JR)*RFAC3*SN(MR+1,INDX)*RFAC4                         02030013
      RCONV3= Q*D(JR)*RFAC3*DSN(MR+1,INDX)*RFAC4                        02040013
      IF (JR.EQ.NR) TN=  RCONV2                                         02050000
      IF (JR.EQ.NR) DTN=RCONV3                                          02060000
      IF (JR.EQ.IFN) TN2=  RCONV2                                       02070000
      IF (JR.EQ.IFN) DTN2=RCONV3                                        02080000
      IF (JR.EQ.IFN) GO TO 8                                            02090000
      R2A=R2A+RCONV2                                                    02100000
      DR2A=DR2A+RCONV3                                                  02110000
 1070 CONTINUE                                                          02120000
    8 IF (MOD(L,2).EQ.0) NR=NR-1                                        02130000
C                                                                       02130100
         CALL  REM(M,NR,H,XI,ALAMD,V,AL,SUM,B,0)                        02140003
         CALL  REM(M,NR,H,XI,ALAMD,DV,ALD,SUMD,BD,1)                    02150003
C                                                                       02150103
      CF=TN/DEXP(-AL*(M+NR))/DFLOAT(M+NR)**B/ DEXP(SUM)                 02160000
      DCF=DTN/DEXP(-ALD*(M+NR))/DFLOAT(M+NR)**BD/ DEXP(SUMD)            02170000
      R2B=(CF*DEXP(-AL*(M+NR+2))/AL*V+TN2)/2.                           02180000
      DR2B=(DCF*DEXP(-ALD*(M+NR+2))/ALD*DV+DTN2)/2.                     02190000
      R2N= XIA*(R2A+R2B)/PI/FCTR                                        02200016
      DR2N= XIB*R2N+C*XIA*(DR2A+DR2B)/PI/FCTR                           02210016
      R2N=R2N/SCLSN                                                     02220000
      DR2N=DR2N/SCLSN                                                   02230000
      DR2NC=(1.0/(C*XI2)+R2N*DR1)/R1                                    02240000
      DFN= DABS((DR2NC-DR2N)/DR2NC)                                     02250000
    7 CONTINUE                                                          02260000
C                                                                       02270003
CC        SPHEROIDAL RADIAL FUNCTIONS OF THE SECOND KIND                02280003
CC        IN TERMS OF THE ASSOCIATED LEGENDRE FUNCTIONS                 02290003
      IF(MOD(L,2).EQ.1) GO TO 3000                                      02300000
      JFIN=2*M                                                          02310000
      JINT=2                                                            02320000
      GO TO 3010                                                        02330000
 3000 JFIN=2*M-1                                                        02340000
      JINT=1                                                            02350000
 3010 QSUM= 0.D0                                                        02360000
      DQSUM= 0.D0                                                       02370000
      DO 10 I=1,IFIN,2                                                  02380000
      IA=I-1                                                            02390000
      IF(MOD(L,2).EQ.1)  IA=I                                           02400000
      MA=M+IA                                                           02410000
      IF(MA.GE.MAXPQ) GO TO 11                                          02420000
      QD=D(I)*DQX(MA+1)                                                 02430000
      QSUM=QSUM+D(I)*QX(MA+1)                                           02440000
      DQSUM=DQSUM+QD                                                    02450000
      IF( DABS(QD/DQSUM).LT.1.0D-13) GO TO 11                           02460000
   10 CONTINUE                                                          02470000
   11 IF(M.EQ.0) GO TO 3030                                             02480000
      DNRAT(JFIN+2)= 0.D0                                               02490000
      DO 20 J=JINT,JFIN,2                                               02500000
      JA=JINT+JFIN-J                                                    02510000
   20 DNRAT(JA)=-E(-JA+2)/(F(-JA)+G(-JA-2)*DNRAT(JA+2))                 02520000
      DN0=D(1)                                                          02530000
      DO 30 JB=JINT,JFIN,2                                              02540000
      DN(JB)=DN0*DNRAT(JB)                                              02550000
      MJB=M-JB                                                          02560000
      MJC=MJB+1                                                         02570000
      IF(MJB.LT.0) GO TO 2050                                           02580000
      QNA=QX(MJB+1)                                                     02590000
      QNB=QX(MJC+1)                                                     02600000
      GO TO 2060                                                        02610000
 2050 MJD=-MJB                                                          02620000
      IF(MJD.EQ.1) GO TO 2070                                           02630000
      QNA=QN(MJD)                                                       02640000
      QNB=QN(MJD-1)                                                     02650000
      GO TO 2060                                                        02660000
 2070 QNA=QN(1)                                                         02670000
      QNB=QX(1)                                                         02680000
 2060 DQN=(DFLOAT(MJB-M+1)*QNB-DFLOAT(MJB+1)*XI*QNA)/XI2                02690000
      QSUM=QSUM+DN(JB)*QNA                                              02700000
      DQSUM=DQSUM+DN(JB)*DQN                                            02710000
      DN0=DN(JB)                                                        02720000
   30 CONTINUE                                                          02730000
      IF(MOD(L,2).EQ.1) GO TO 3040                                      02740000
      COEF2=FF(M,N)*FACTMM(M)*DN(JFIN)*C*PI/(DFLOAT(M+M-1)*C**M)        02750013
     +      *FCTR                                                       02760013
      GO TO 3050                                                        02770000
 3040 COEF2=-FF(M,N)*FACTMM(M)*DN(JFIN)*C2*PI/(DFLOAT((M+M-3)*(M+M-1))* 02780013
     +       C**M)*FCTR                                                 02790013
 3050 CONTINUE                                                          02800000
      GO TO 3060                                                        02810000
 3030 IF(MOD(L,2).EQ.1) GO TO 3070                                      02820000
      COEF2=-FF(M,N)*D(1)*C*PI                                          02830000
      GO TO 3060                                                        02840000
 3070 COEF2=-FF(M,N)*D(1)*C2*PI/3.0D0                                   02850000
 3060 CONTINUE                                                          02860000
      IH=H                                                              02870000
      IF(MOD(L,2).EQ.1) GO TO 3080                                      02880000
      KINT=M+M+2                                                        02890000
      KFIN=30+2*(M+IH)                                                  02900000
      GO TO 3090                                                        02910000
 3080 KINT=M+M+1                                                        02920000
      KFIN=31+2*(M+IH)                                                  02930000
 3090 IF(KFIN.GE.138) KFIN=138+MOD(L,2)                                 02940011
 2040 DPRAT(KFIN+2)= 0.D0                                               02950000
      DP0=D(1)                                                          02960000
      IF(M.NE.0)  DP0=DN(JFIN)                                          02970000
      DO 40 KA=KINT,KFIN,2                                              02980000
      KB=KINT+KFIN-KA                                                   02990000
      IF(KB.EQ.KINT) GO TO 2010                                         03000000
      DPRAT(KB)=-E(-KB+2)/(F(-KB)+G(-KB-2)*DPRAT(KB+2))                 03010000
      GO TO 40                                                          03020000
 2010 IF(MOD(L,2).EQ.1) GO TO 2020                                      03030000
      DPRAT(KB)=C2/FLOAT((M+M+1)*(M+M-1))/(F(-KB)+G(-KB-2)*DPRAT(KB+2)) 03040000
      GO TO 40                                                          03050000
 2020 DPRAT(KB)=-C2/FLOAT((M+M-1)*(M+M-3))/(F(-KB)+G(-KB-2)*DPRAT(KB+2))03060000
   40 CONTINUE                                                          03070000
      PSUM= 0.D0                                                        03080000
      DPSUM= 0.D0                                                       03090000
      DO 50 KC=KINT,KFIN,2                                              03100000
      DP(KC)=DP0*DPRAT(KC)                                              03110000
      KD=KC-M-1                                                         03120000
      KE=KD+1-M                                                         03130000
      IF(KE.GT.MAXPQ) GO TO 50                                          03140000
      PSUM=PSUM+DP(KC)*PX(KE)                                           03150000
      DPSUM=DPSUM+DP(KC)*DPX(KE)                                        03160000
      DP0=DP(KC)                                                        03170000
   50 CONTINUE                                                          03180000
      KFINM=KFIN-M-M                                                    03190000
      IF(KFINM.GE.MAXPQ) GO TO 2030                                     03200000
      IF(KFIN.GE.137) GO TO 2030                                        03210011
      CRIT= DABS(DP(KFIN)*PX(KFINM)/PSUM)                               03220000
      KFIN=KFIN+4                                                       03230000
      IF(CRIT.GT.PSCRT) GO TO 2040                                      03240000
 2030 KFIN=KFIN-4                                                       03250000
      R2Q=(QSUM+PSUM)/COEF2                                             03260000
      DR2Q=(DQSUM+DPSUM)/COEF2                                          03270000
      DR2QC=(1.0/(C*XI2)+R2Q*DR1)/R1                                    03280000
      DFQ= DABS((DR2QC-DR2Q)/DR2QC)                                     03290000
      IF(XI.LT.1.005) GO TO 1                                           03300000
      IF(HXX.LT.10.0) GO TO 1                                           03310000
      IF(L.GT.LMTN) GO TO 1                                             03320000
      IF(DFQ.LE.DFN) GO TO 1                                            03330000
      R2=R2N                                                            03340000
      DR2=DR2N                                                          03350000
      GO TO 4                                                           03360000
    1 R2=R2Q                                                            03370000
      DR2=DR2Q                                                          03380000
      GO TO 4                                                           03390000
    3 WRITE(6,100)                                                      03400000
  100 FORMAT(1H0,20X,'*  RPRSWF  -  RETURN WITHOUT CALCULATION  *  COEFF03410000
     1ICIENTS D ARE ALL ZERO' )                                         03420000
    4 RETURN                                                            03430000
      END                                                               03440000
